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The spectrum of baryons containing three b quarks is calculated in nonperturbative QCD, using 
the lattice regularization. The energies of ten excited bbb states with = > 
and I are determined with high precision. A domain- wall action is used for the up-, down- and 
strange quarks, and the bottom quarks are implemented with NRQCD. The computations are done 
at lattice spacings of a « 0.11 fm and a ~ 0.08 fm, and the results demonstrate the improvement of 
rotational symmetry as a is reduced. A large lattice volume of (2.7fm)"^ is used, and extrapolations of 
the bbb spectrum to realistic values of the light sea-quark masses are performed. All spin-dependent 
energy splittings are resolved with total uncertainties of order 1 MeV, and the dependence of these 
splittings on the couplings in the NRQCD action is analyzed. 

PACS numbers: 12.38. Gc, 14.20.Mr 



I. INTRODUCTION 



Heavy quarkonium has been studied in great detail both experimentaUy and theoretically. Because its valence 
quark masses are large compared to Aqcd, heavy quarkonium is an excellent system for probing the strong force on 
multiple scales [T]. In addition to these familiar heavy quark-antiquark bound states, QCD also predicts the existence 
of an analogous system in the baryonic sector: the bound states of three heavy quarks. Given the huge importance 
of quarkonium, it is desirable to investigate triply-heavy baryons in similar depth. 

Most of the existing continuum-based calculations of triply-heavy baryon masses used rather simple models [2H20] . 
Correspondingly, the results have uncontrolled uncertainties and are widely spread. A more systematic continuum- 
based approach is potential NRQCD (pNRQCD), which was formulated for baryons containing two or three heavy 
quarks in Ref. ^I]. In combination with the NNLO perturbative static three-quark potential [22], pNRQCD was 
recently used to estimate the ccc, ccb, ebb, and bbb masses [23| . 

No experimental results are available so far for triply heavy baryons (see Ref. 24J for a recent calculation of 
production cross sections at the LHC). This means that first-principles nonperturbative lattice QCD calculations are 
essential to test the model-dependent approaches. For the ^bbb, the ground-state mass was already calculated with 
high precision using lattice QCD in Ref. |2S]. However, much more information about the interactions between three 
heavy quarks can be gained by computing the spectrum of the corresponding excited states, including in particular 
the spin-dependent energy splittings. The first such calculation of bbb excited states in lattice QCD is reported here. 
Lattice calculations of light-baryon excited states can be found for example in Refs. p6l434j . 

To fully accommodate the physics of the light sea quarks in lattice QCD, the spatial box size L has to be chosen 
such that L ^ I/tott. With the presently available computing resources, this requirement means that the lattice 
spacing is too coarse to treat the b quarks in the same way as the light quarks. Therefore, as in Ref. [5S], the b quarks 
are implemented here with improved lattice NRQCD [351 [35]. NRQCD is an effective field theory for heavy quarks 
that retains all the gluon and light-quark degrees of freedom without change. For the heavy quark Lagrangian, a 
nonrelativistic expansion is performed in powers of the heavy-quark velocity v, and the coefficients of the NRQCD 
effective operators are determined by matching to QCD. Thereby, the results of QCD can be reproduced in principle 
to an arbitrary order in v. For bb and bbb hadrons, one has (u^) ~ 0.1. The lattice NRQCD action used in Ref. [2S] 
was complete through order u^. Because the present work aims to accurately compute also spin-dependent bbb energy 
splittings (fine and hyperfine structure), here the spin-dependent order-u^ terms are included in the NRQCD action, 
as already done in the calculation of the bottomonium spectrum of Ref. [37]. Furthermore, the coefficients of the 
leading spin-dependent operators, which are of order w^, are tuned nonperturbatively. 

As usual in lattice QCD, the Euclidean path integral is performed by averaging over importance-sampled gauge 
field configurations. The ensembles of gauge fields used here match those used in Refs. [3S] and [37], and have been 
generated by the RBC/UKQCD collaboration [3S]. These ensembles include the effects of dynamical u-, d- and 
s- quarks, which were implemented using a domain- wall action i39ri4T] . Seven different ensembles with a range of 
light-quark masses and lattice spacings of a w 0.11 fm and a ss 0.08 fm are included in the analysis. 

The bbb energy levels are extracted from the time-dependence of Euclidean two-point functions of interpolating 
operators with the desired quantum numbers. The construction of these interpolating operators, which takes into 
account the reduction of the continuum rotational symmetries to the lattice rotational symmetries, follows the highly 
successful method originally developed for light baryons in Ref. '31'. This method, as well as the computation of the 
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bbb two-point functions, is explained in Sec.|TTj The details of the lattice actions and parameters are given in Sec. Ill 



Next, Sec. IV describes the fitting of the two-point functions and the angular momentum identification. The spectrum 
results are extrapolated in the light-quark masses to obtain the final results in Sec.|Vj An additional section (Sec. VI) 
is devoted to investigating the dependence of the bbb energy splittings on the various operators in the NRQCD action. 



II. CONSTRUCTION OF BARYON INTERPOLATING OPERATORS 



In this section we construct interpolating operators, fl, that give access to bbb states up to J = |. The method is 
taken from Ref. [34] j but is described again in the following specifically for the case needed here, where all three quark 
fiavors are equal and only two-component Pauli spinors are used. Going through the derivation of the interpolating 
operators also gives some insight into the structure of the bbb states extracted in the numerical part of the calculation. 
However, it is important to remember that the spectrum calculated here is that of the (lattice) QCD+NRQCD 
Hamiltonian: Hqcj^ln) = En\n). The interpolating operators determine only the overlap factors (n|r2|0), not the 
energies E„. For the numerical calculation it is nevertheless advantageous to construct operators that have large 
overlaps only with selected bbb states, to get good statistical precision for the energy levels and identify their angular 
momentum J. 

A key feature of the approach from Ref. |34j is the initial construction of operators with definite quantum numbers 
J and m according to the continuum rotational symmetry (Sec. II A). This is then followed by the subduction, where 
linear combinations of the different m-components at a given J are formed such that these transform irreducibly under 
the lattice rotational symmetries (Sec. II B). The numerical calculations demonstrate that the rotational symmetry 
breaking is very weak, and operators subduced from continuum operators with different values of J retain an ap- 
proximate orthogonality even if they fall in the same irreducible representation of the octahedral group. This feature 
dramatically simplifies the angular momentum identification for the extracted energy levels. 

Following the group-theoretical operator construction. Sec. II C then describes the initial smearing of the quark 
fields and the calculation of the baryon two-point functions on the lattice. 



A. Operators with definite continuum J 



In all baryon operators, the colors of the three quarks are combined into a singlet using the totally antisymmetric 
color wave function Cahc- In the case considered here, the three quarks have equal fiavor. Therefore, to satisfy the 
Pauli principle, the product of the spin and spatial wave functions must be totally symmetric. The spatial structure is 
obtained by applying up to two derivative operators to Gaussian-smeared quark fields. The derivatives are combined 
to a definite total orbital angular momentum L and a definite permutation symmetry. Similarly, the spins of the 
three quarks are combined to a definite total spin S and definite permutation symmetry. Finally the derivative and 
spin wave functions obtained in these two separate steps are combined to obtain baryon operators with a definite 
total angular momentum J and the desired total symmetry of the product of the spin and spatial wave functions. 
Note that L and S are not conserved quantum numbers, and are only used to label the structure of the interpolating 
operators. 

We begin by combining the three quark fields to definite total spin S. Because NRQCD is used for the heavy 
quarks, there are only two spin components, denoted by ip^ and ip^. The color indices are omitted here, but remain 
uncontracted at this stage (the contraction with tabc is only performed after the gauge-covariant derivatives have been 
applied). The »S' = | combinations are given by 
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where the subscript 5 indicates the total symmetry under permutations. For S* = ^, one can construct both mixed- 
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symmetric (MS) or mixed-antisymmetric (MA) combinations: 



Oma(|, -5) -;^(V';'0tV'; - V'tV';'^;)- (3) 

Next, we come to the derivatives. A single derivative is an L = 1 object, with m-components given by 

D±i=±\{D^±iDy), 
D.^-^,D^. (4) 



Recall that in this section we work in continuous space; lattice derivatives will be defined in Sec. |II C| In the following, 
we use the notation Dm^ for a derivative acting on the fc-th quark in the baryon operator. As in Ref. [34j, we define 
the following combinations with definite permutation symmetry: 

i?W(l,m) = -^(2i?(3)-i?W-i?(„2)), 

i^W (i.H = ^(^ii^-i^L'')- (5) 



(No totally antisymmetric combination exists, and the totally symmetric combination vanishes at zero-momentum.) 
Using the Clebsch-Gordan coefficients (L, to|1, mi; 1, TO2), we can combine two single-derivative operators of the form 
([5]) into double-derivative operators with definite total L and definite permutation symmetry as follows [34] : 
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-D[,^l,(l,mi)DW(l,m2)) 


Df{l,m) 




l,TOi; l,m2)(^ 


+ i?W(l'"^l)^Mk(l-"^2)- 


i5[Jl,(l,mi)i?W(l,m2)). 



(6) 

The first three of the above combinations can give either L = or L = 2, while the last combination is restricted to 
L = 1. 

Now we can combine the spin- and spatial wave functions, distinguishing the cases of zero, one, and two derivatives. 
Without derivatives, the requirement of total symmetry restricts the spin to 5 = |. Since L = 0, we only get J — 
in this case: 

[OsiWf'^Osilm). (7) 

In one-derivative baryon operators, the derivative part, Eq. ([5]), always has mixed symmetry. Therefore, to get 
a totally symmetric combination, the spin part must also have mixed symmetry, and hence S = ^. Because the 
derivative has L = 1, we can combine L and 5* to the total angular momenta J — \ and J = |: 

[Z)W(l) OM(i)];l=^'^ = ^ (•^>™|l:'™i;5."^2)(i?^L(l'"^i)OMs(im2)-t-i^^Ul'"^i)OMA(i,m2)). (8) 



mi ,m2 



Finally, we consider the double-derivative operators. Because no totally antisymmetric spin combinations exist, the 
totally antisymmetric derivative combination in the last line of Eq. ([6]) is excluded, and the two derivatives can only 
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combine to L = or L = 2. In each case, one can have S = with mixed symmetry or = | with total symmetry. 
Thus, one obtains the following combinations: 



[Df\0) Os(i)];-2 =4'k0,0)Os(i,m) 



,/=; 



3 5 



13 5 7 



—Jk X! (-^i '7i|2,mi; i,m2)(^L'Ms(2j'™i) Oms(5,?772) +i?MA(2,'7ii) Oma(5,"12)^ 

^ mi, 7712 

^ (J,m|2,mi;|,TO2)4''(2,mi)Os(i,rn2). (9) 



f2l 

Note that the combination with (0, 0), which corresponds to the spatial Laplacian, was excluded in Ref 
the argument that it vanishes at zero momentum. However, this is no t the case for the method of smearing the quark 

r [21 „ , r- 

fields and constructing the two-point functions described in Sec. II C In fact, the operator [Dg (0) 
good overlap with the first radially excited J = \ state, and including this operator in the basis significantly improves 
the extraction of this energy level. 



with 
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B. Subduction to irreducible representations of the double cover of the octahedral group 

In the previous section, we constructed operators \p\ ^ that transform under rotations like the basis vectors | J, m) 
of irreducible representations of SU{2). The group SU{2) is the double cover of the continuum three-dimensional 
rotation group S0{3). On a cubic lattice, the rotational symmetry is reduced to the discrete group ^O, the double 
cover of the octahedral group O. The group ^O, which is obtained from O by adding a negative identity for ±2tt 
rotations, has 48 elements in 8 conjugacy classes. Correspondingly, has 8 irreducible representations denoted as 
Ai, Ai, E, Ti, T2, Gi, G2, H (see, for example, Ref. |42J)- Their dimensions are 1, 1, 2, 3, 3, 2, 2, 4, respectively. 

Starting from an operator [f^],^, it is possible to form suitable linear combinations of its different m-components, so 
that these linear combinations transform in irreducible representations. A, of the double-cover octahedral group: 

m 

This process is referred to as reduction or subduction [341 142] . and the coefficients form the subduction matrices. 
Here, "A denotes the n-th occurrence of an irrep A of ^O, and r = 1, dim(A) denotes its row (like m denotes the 
row for the SU{2) irrep). For each value of J, only selected irreps of appear in the subduction, such that the 
sum of their dimensions equals 2 J -I- 1 (the dimension of the original SU(2) irrep J). This is indicated in Table |I] 
For integer values of J, only the irreps Ai, A2, E, Ti, and T2 appear. Conversely, for half-integer J, only the irreps 
Gi, G2, and H occur. Since we are considering baryons, we will only be concerned with these three irreps in the 
remainder of the paper. The subduction matrices for ( J = — > Gi and ( J = §) — > iJ are simply the 2x2 and 4x4 

identity matrices, so that, for example, \P\q ^ = ^ = The subduction matrices for J = | 

and J — \ can be found in Ref. |34j . 

So far we have only discussed the rotational symmetry. Additionally, we can classify the operators according to 
their transformation properties under space inversion, which remains an exact symmetry on the lattice. Then all of 
the irreducible representations come in parity-even and parity-odd versions, as indicated by subscripts g (gerade) and 
u (ungerade): Aig, T2g, Gig, G2g, Hg, and T2U, G2t,, H^- In this work, the baryon operators are 

constructed from two-component NRQCD spinors, and therefore the parity of an operator is determined entirely by 
the number of derivatives it contains: an even number of derivatives corresponds to even parity and an odd number 
of derivatives corresponds to odd parity. 

The 11 different baryon operators constructed in Eqs. 
each in the Gig and G2g irreps, and 1 operator each in the Gi^ and irreps. This set of operators is summarized 
in Table HD 
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TABLE I. Subduction of SU{2) irreps to irreps, up to J = | (from Ref. [42]). 
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TABLE IL Interpolating operators, named according to their parity (cj: +, u: — ) and irreducible representation of ^O. The 
superscript labels the different operators within a given irrep and parity. 



C. Computation of two-point functions on the lattice 



The group-theoretical construction of baryon operators through subduction was performed here in the same way as 
done for light baryons in Ref. [34]. However, the method for smearing the quark fields and computing the two-point 
functions in terms of quark propagators differs from that used in Ref. [34] . Instead of distillation |43j , here the more 
traditional approach starting from Gaussian-smeared point sources, as in Ref. [H], is chosen. This has the advantage 
over distillation that the quark smearing width can be made very narrow without increasing the computational cost. 
A narrow smearing width is needed to get a good overlap with the physical bbb states, which are expected to be very 
small objects as a consequence of the large &-quark mass. 

In the approach used here, the smeared 6-quark fields ^ entering in Eqs. 
quark fields ip through 

/ 7-2 \ 



ip) 



are defined in terms of the unsmeared 
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where A^^^ is a three-dimensional gauge-covariant lattice Laplace operator, 

1 ^ 

A(2)i/,(a;,i) = -—J2[Uj{x,t)i;{x + aj,t)^2^P{x,t) + U^j{x,t)^{x^aj,t) 



(12) 



In this work, a smearing radius of rg « 0.14 fm is used in Eq. (111. The gauge-covariant derivatives in the baryon 
operators then act on these smeared quark fields. The continuous derivatives Dj used in Sec. II A are replaced by 
lattice versions Vj, which are defined as 



2a 



Uj{x, t)ip{x + aj, t) — U-j{x, t)ip{x — aj, t) 



(13) 



The tilde on the gauge links in Eqs. ( 12 ) and ( 13 ) indicates that these are also smeared, using the procedure of 
Ref. [45 . The gauge link smearing in the hadron interpolating fields is performed to reduce statistical noise |44| . 
The baryon operators constructed in the previous two sections contain quark fields with up to two derivatives. It is 
convenient to introduce new objects ipi, where i labels all the required thirteen derivative combinations: 
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(14) 

Additionally to the derivative index i — 1, ...13, these fields -0^ — (ipaai) have a color index a — 1,2,3 and a spinor 
index a — 1,2 {— ],). Then, all baryon interpolating operators used here have the form 



^rix, t) =Taipjjk <^abc '4'aai{x,t) ■ll!hi3jix,t) tl>cjk{x,t), 



(15) 



where Tai is the set of complex- valued coefficients from Sec. II B for each operator. An example is given in Table 



III The two-point function at zero momentum, allowing different operators VIy and fir' at sink and source, is then 
defined as 



CT,v'{t t') = ^{nr{x,t) nl,{x',t') 

X 

= '^^^aiPj-ik ^abc ^aiPjjk^abc 



X {i^aai{x,t) i;bpj{x,t) ^c^k{x,t) ^l_.{x',t') -j/ii^. (oj' , t') 0^_^(a;',t') 



(16) 



where the brackets denote the Euclidean path integral over the gauge fields and fermions, weig hted with e"-^. The 
path integral over the fermions can be performed explicitly, giving heavy-quark propagators and determinants of the 
Dirac operators for all quark flavors. Following Ref. ^], we define three-quark propagators (for a given gauge field 
U) that have been color-contracted and summed over x: 



G 
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ai ai Pj (3j jk jk 



, ;^j,(^j ^ J 3; ) — Eabc ^abc Gaaiaai{x, t,X ,t) Gj^pj i,pj{x,t,X ,t ) G ^^p. c^k{x, t,X ,t), (17) 



where Gg^„ig^g,j{x,t,x' ,t') denotes a heavy-quark propagator with smearing and, depending on i and i, derivatives 
at the source and sink. Performing the fermionic path integral in Eq. ( 16 1 gives six contractions because all three 
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TABLE III. Example for the baryon operator coefficients Taiisj-yk- Shown here are the non-zero coefficients for the first row 
(r = 1) of the operator Hi^\ 



heavy-quark flavors are equal. Using the antisymmetry of the epsilon tensor, one obtains 

Cr, r' (i - t') Tai pj jk T^*. ^aLl « 0] -yk 7fc ^' 

~^ ^ ai jk ai ^ ' 
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(18) 



Here, ( ... )u denotes the path integral over the gauge fields U only, where the weighting factor is given by e"'^'!""*" x 
(fermion determinants). 

In the numerical calculations, performing all the multiplications in the three-quark propagator (171 is expensive, 
and it is important to use symmetries to reduce the number of operations needed. Defining multi- indices I = {ai ai), 
J = [jSj f3j), and K = {"fk^k), one finds that G^^jj^ is totally symmetric in /, J, K. Furthermore, since the baryon 
operators constructed in the previous two sections contain at most two derivatives total, only those components of 
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with 



noii) + noU) + noik) < 2, ?in(i) + noij) + noik) < 2 



(19) 



are needed [n]j(i) denotes the number of derivatives associated with the index i, see Eq. (14)]. 



III. LATTICE ACTIONS AND PARAMETERS 



The path integral over the gauge fields U in Eq. (18) is performed by averaging over samples of lattice gauge field 
configurations. The configurations used here have been generated by the RBC/UKQCD collaboration [3S], and include 
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0.8439 


1.196 


1.168 


1480 - 


8520, 10 


32 


0.1139(19) 


0.4194(70) 


24^ 


X 64 


2.13 


0.02 


0.04 


2.622 


0.8433 


1.196 


1.168 


1800 - 


3600 10 


32 


1177f29'l 


541('14~) 


24^ 


X 64 


2.13 


0.03 


0.04 


2.691 


0.8428 


1.196 


1.168 


1280 - 


3060, 10 


32 


0.1196(29) 


0.641(15) 


32^ 


X 64 


2.25 


0.004 


0.03 


1.831 


0.8609 


1.175 


1.113 


580 - 


6840, 10 


24 


0.0849(12) 


0.2950(40) 


32^ 


X 64 


2.25 


0.006 


0.03 


1.829 


0.8608 


1.175 


1.113 


552 - 


7632, 16 


24 


0.0848(17) 


0.3529(69) 


32^ 
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TABLE IV. Summary of lattice parameters. The coupling in the Iwasaki gauge action is given as /3 — d/g^ , and amu,d, arUs 
are the bare masses of the domain-wall sea quarks. The parameters ami,, uql, cs, and C4 enter in the NRQCD action for the 
b quarks. The lattice spacings, a, were computed in Ref. [37]. The MD (molecular dynamics) range specifies the range of the 
gauge-field generation Markov chain [HH] for which "measurements" are performed. Th e m easurements are separated by the 
given step size in MD time, and are done for Usic different source locations [(a;', t') in Eq. ( 16 1] on each gauge field configuration. 



dynamical u-, d-, and s- quarks, with m„ = m^. These quarks were implemented with a domain-wall action |39l - 
|4T] , which is a five-dimensional action that leads to an approximate lattice chiral symmetry for the four-dimensional 
theory. This chiral symmetry becomes exact when the length of the auxiliary fifth dimension is taken to infinity. For 
the gauge action, the Iwasaki discretization [JSIHT] is used (the gauge fields are four-dimensional, i.e. constant in the 
5-direction). The domain- wall formalism requires additional Pauli-Villars fields to cancel bulk modes [101 HH], so that 
the gauge fields U are distributed with probability density proportional to 



det[KPW(^. aM5,amu,d)? det[K^'^{U; aM5,am,)] _s,,„,,[c/] 
det[K^^iU;aM5,lW e ^ ^ , 



(20) 



where K^^{U]aM^ , am) is the five-dimensional domain- wall operator with domain-wall height M5 and quark-mass 
TO. Seven ensembles of gauge fields with different parameters are included in the analysis, as shown in Table |IV| 
There are ensembles with two different values of the gauge coupling /3 = 6/(7^, leading to lattice spacings of a ss 0.11 
fm and a ~ 0.085 fm. The number of lattice points is chosen to be 24^ x 64 and 32'^ x 64, respectively, so that the 
spatial volume in physical units is equal to about (2.7 fm)"^ in both cases. 

The lattice NRQCD action for the h quarks has the same form as in Ref. [37j. It can be written as 

Stp = a 

t) t) - K{t) ij{x, t-a)], (21) 
where t/j is a two-component spinor, and K(t) is given by |36] 



with the leading-order kinetic energy operator. 
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and the following higher-order relativistic and discretization corrections: 
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Here, E and B are the chromoelectric and chromomagnetic components of a lattice gluon field strength tensor. Unlike 
in the previous sections, the tilde appearing on some of the quantities in Eq. (24) does not denote smearing; instead 
it denotes improvement corrections which reduce discretization errors [33]. The action is also tadpole-improved |49j . 
with the values of the Landau gauge mean link uql as given in Table |IV[ The heavy-quark masses in lattice units, 
anih, are set to the physical values as determined for the same gauge field ensembles in Ref. [37j. 

In Eq. (24), the terms with matching coefficients ci, C2, C3, and C4 are the relativistic corrections of order v'^. The 



terms with coefficients C5 and cg are spatial and temporal discretization improvements for Hq. Finally, the terms with 
coefficients cy, cs, and cg are the spin-dependent order-u® terms. In principle, additional operators containing four 
(or more) quark fields are introduced through gluon loops, but these are not included here. 

At tree level in the matching of NRQCD to QCD, the coefficients Ci in Eq. (24) are all equal to 1 
in SH are suppressed relative to Hq by at least one power of 

accuracy of order agv'^ « 0.02 for the radial and orbital energy splittings in the bb and bbb systems 



Because the terms 
using the tree-level values for q already provides 

However, spin 



splittings first arise through the operators with coefficients C3 and C4, and therefore these two coefficients are tuned 
nonperturbatively here. The tuning condition used here is that, when calculated with the lattice NRQCD action, the 
following two combinations of bottomonium IP energy levels agree with experiment: 



(25) 



and 



mxbo) + mxbi)-E{xb2). 



(26) 



As discussed in Ref. |50| and confirmed numerically in Ref. [37j, to a good approximation the combi nat ion ( 25t i s 
proportional to C3, while (26) is proportional to c|. Table VII of Ref. [57] gives numerical results for (25) and (|26|), 
computed with Ci — 1 for the same order-u^ NRQCD action on the same gauge field ensembles. Using these results, 
one can then solve for C3 and C4 so that the experimental values [5T] for (25) and (26) are reproduced: 



C3 



1.196 ± 0.106, a « 0.11 fm, 
1.175 ±0.084, a « 0.08 fm, 



C4 



1.168 ±0.081, a; 
1.113 ± 0.053, a; 



0.11 fm, 
0.08 fm. 



(27) 



In the present work, the main calculations of the bbb spectrum are performed directly at C3 and C4 set equal to the 
central values in Eq. (27 1, and with ci — C2 = ~ cq ~ — cs — cg — 1. The uncertainties in (|27[) are mainly 



statistical, and the resulting uncertainties in the bbb spectrum will be included in the final results (Sec. |V[) 



IV. FITS OF THE TWO-POINT FUNCTIONS AND ANGULAR MOMENTUM IDENTIFICATION 



The two-point functions defined in Eq. (16) are labeled by F and F', which determine the baryon interpolating 
operators at the sink and source, respectively. The two-point functions vanish when F and F' correspond to different 
irreducible representations (irreps) of the double-cover octahedral group, or when F and F' correspond to different 
rows of the same irrep. In the remaining cases of equal irrep and equal row at source and sink, one can average over 
the different rows. In the following, we use the notation A^*^ for row r of the i-th operator in irrep A, according to 
Table [H] Then the row-averaged two-point functions are defined as 



dim(A) 

^^(^-0-^ E ^A(^AO)(^-*')- (28) 

For the operators in Table |llj one obtains a (7 x 7) matrix of two-point functions in the Hg irrep, (3 x 3) matrices 
in the Gig and G2g irreps, and (1 x 1) "matrices" in the Hu and Giu irreps. The magnitudes of the rescaled two- 
point functions \c[y'\/ \J~c\^^~C^ at one time slice are shown in Figs.llj j2j andjsjfor the Hg, Gig, and G2g irreps, 

between ope: 



respectively. The first important observation is that cross-correlations between operators subduced from continuum 
operators that differ in at least one of the quantum numbers L, 5*, or J are small. Note that J is an exactly conserved 
quantum number in the continuum, but L and S are not. The weak coupling between operators subduced from 
different J-values indicates that rotational symmetry breaking by the lattice is small. This has also been observed 
in Ref. ^34] for light baryons. On the other hand, the weak coupling between operators subduced from common 
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FIG. 1. Visualization of rescaled matrix two-point functions \Cij\/ y/CuCjj in the Hg irreducible representation, at one time 
slice. Off-diagonal entries larger than 0.01 are also given numerically (the i = 1, j = 2 entry is 0.98). Left plot: a ~ 0.11 fm, 
am^^d = 0.005, {t - t')/a = 5. Right plot: a ^ 0.08 fm, am„,d = 0.004, (t - t')/a = 6. 
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FIG. 2. Visualization of rescaled matrix two-point functions \Cij\/ ^CaCjj in the Gig irreducible representation, at one time 
slice. Off-diagonal entries larger than 0.01 are also given numerically. Left plot: a « 0.11 fm, am^.d ~ 0.005, {t — t')/a = 5. 
Right plot: a ^ 0.08 fm, am„,d = 0.004, {t - t')/a = 6. 
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FIG. 3. Visualization of rescaled matrix two-point functions \Cij\l ^CaCjj in the G2g irreducible representation, at one time 
slice. Off-diagonal entries larger than 0.01 are also given numerically. Left plot: a ~ 0.11 fm, amu,d ~ 0.005, (t ~ t')/a = 5. 
Right plot: a ^ 0.08 fm, am„,d = 0.004, {t - t')/a = 6. 



J-values but different L- or different 5-values is a new feature appearing here. Because of the large mass of the 
b quarks, the dynamics is approximately nonrelativistic, and the spin-orbit coupling is suppressed, so that L and 
S are approximately conserved. In fact, for the lattice spacings considered here, the operator overlaps between 
different L- or S'-values appear to be smaller than that between different J-values. Furthermore, the overlaps between 
operators subduced from different J-values (for example between Hg^^ and Hg^"^ , which are subduced from J = | and 
J = |, respectively) appear to be somewhat larger than what was seen for light baryons in Ref. [34] . This may be 
a consequence of the much smaller physical extent of the bbb baryons [as modelled by the initial smearing width of 
rs ~ 0.14 fm in Eq. ([lT|], which makes the operators more sensitive to the non-zero lattice spacing. 

As can be seen in Fig. [l] there is a strong overlap between the Hg^^ and iJg^-* operators, because both are subduced 
from continuum operators with the common quantum numbers L = 0, S = |, J = |. All other cross-correlations, 
also in the Gig and irreps (Figs. [2] and |3| are small, because there is suppression as a consequence of different J, 
L, or S. 
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Further information can be gained by looking at the lattice-spacing dependence of the operator overlaps. In each 
of the figures, the left plot shows data from a « 0.11 fm, while the right plot shows data from a ~ 0.08 fm. It can be 
seen that the cross-correlations between operators subduced from different continuum J are smaller at the finer lattice 
spacing, demonstrating the improvement of rotational symmetry as a is reduced. On the other hand, the overlaps 

(3) fl) (3) (2) 

between Hg and Hg , as well as between Hg and Hg , are not smaller at the finer lattice spacing. In that case, 
the operators are all subduced from the same J {— |), and one does not expect the cross-correlations to vanish in 
the continuum limit. 

In this work, the matrix two-point functions in each irrep A were fitted directly using the form 

N 

C^it-i!) = Y.<Kl^''^'^-''- (29) 

The number of exponentials was chosen to be equal to the dimension of the matrix, i.e. equal to the number of 
interpolating operators for each irrep: TV = 7 for i/g, = 3 for Gig and G2g, and = 1 for and Giu- Of course, 
the complete spectral decomposition of the two-point functions also contains an infinite number of higher-energy 
exponentials. Therefore, only the data with t — t' > t^i^ with sufficiently large t^i^ were included in the fit, so that 
the contributions from these higher states are negligible. The dependence of the results on tmin will be discussed later. 

The fits performed here fully take into account the statistical correlations between all data points. The dimension 
of the data correlation matrix for an (TV x iV) matrix fit is equal to Nt N'^, where Nt is the number of time slices 
included in the fit (Nt — imax/a~imin/a + l)- The definition of contains the inverse of this data correlation matrix, 
and one has to make sure that the number of measurements used to estimate the data correlation matrix is much 
larger than its dimension. Because the number of measurements was of order rigrc x ndg ~ 10^ for each ensemble, 
these large, fully correlated matrix fits were possible here (for sufficiently small Nt). In order to reduce the dimension 
of the data correlation matrix to Nt N{N + l)/2 and thereby allow slightly larger Nt, the symmetry of the data in i, j 
(which is exact for infinite statistics) was used. The data for the two-point functions were first symmetrized explicitly 



measurement by measurement, and then the fits using Eq. (29) were performed only for i > j. 

Within each irrep A, the operators A^'^ in Table \n\ are labeled by i such that they are ordered by the energy of the 
state with which they have the strongest overlap (this ordering was not known a priori and was only assigned after 



some initial fits). For each irrep A, the amplitudes in Eq. (29) were then rewritten as follows 



4(A) _ j for n = i, 

forn^z, ^''^ 

using the new parameters v4,[^^ and B^f^^ instead of A^^^' in the fits. The parameters B^l^^ then describe the overlaps 
of the operator A^*-' with the other states n ^ i, relative to the state with n = i. 



Furthermore, the energies -EI'^'' in Eq. i29l were rewritten for n > 1 as 



Ei^^^E^^^+6[^'U...+Si% with (31) 

using the ground-state energy e[^^ and the energy splittings 5[^\ (jjy^lj^ (all in units of 1/a) as the actual fit 
parameters. When computing E^'' (and other combinations of energy levels) from the fit results for e[^^ and S[^\ 
•■■I ^Tv-^-ij the uncertainties were added in a fully covariant way, using the parameter covariance matrix obtained from 
the second derivatives of x^- 

Following Ref . [34] , the spectral overlaps A^^J are used here to assign values of the continuum angular momentum 

J to each energy level E'i'^'*. Examples of fitted energies together with the relative overlap factors A^j^J/aI'^\ 

are shown in Fig.[4]for the Hg, Gig, and G2g irreps (in the cases of the Gi„ and irreps, there is only one operator 
each, subduced trivially from J — and J = |, respectively). The angular momentum identification proceeds as 

follows: for each energy level Ei^\ the operator A^*) with the largest relative overlap factor A^^j/Af^'' is determined. 
The value of J from which this operator was subduced is then assigned to this energy level. As can be seen in Fig. [4j 
no ambiguity arises here. Notice that the two J — ^ levels appearing in the Hg irrep also show up in the G2g irrep, 
with nearly identical energies. Similarly, the J — level appears in all three irreps Hg, Gig, and G2g, again with 
nearly identical energies. For these levels, the absolute overlap factors were also found to be consistent across the 
different irreps, confirming the assignment of J. 

Because of the strong statistical correlations across irreps, the tiny splittings of the J = ^ and J = ^ levels into the 
different lattice irreps, which are caused by rotational symmetry breaking, can be computed with smaller uncertainties 
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FIG. 4. Fitted energies E^f^^ (in lattice unit s; f rom b otto m to top: n = 1, A'^), together with histograms of the corresponding 
relative overlap factors A^^^/A^^^ [see Eqs. (29 1 and (30l]. The fits for the three different irreps were performed independently. 



For each i, the continuum angular momentum J from which the operator A'*' was subduced is given at the bottom. These 
values of J are also indicated by the colors in the histograms (red: J — ^, green: J = |, blue: J = |, orange: J = |). The 
data shown here are from the ensemble with a ~ 0.08 fm and amu,d ~ 0.004; the fits have Wn/a ~ 6. 



than the individual energies of these levels. To this end, simultaneous fits of the two-point functions in the Hg, Gig 
and C?2g irreps were performed, where a global correlated was formed but all fit parameters remained independent 
for each irrep. The results for the rotational-symmetry-breaking-induced energy splittings, converted to MeV, are 
given in Table |V] for two gauge field ensembles. Up to some statistical fluctuations, the splittings are smaller at 
a 0.08 fm compared to a 0.11 fm, consistent with the discretization errors proportional to that are expected 
for the improved lattice NRQCD action used here. Along with the behavior of the off-diagonal matrix elements that 
was discussed at the beginning of this section, the results shown in Table |V] provide another demonstration of the 
improvement of rotational symmetry when the lattice spacing a is reduced. 

Finally, to get the best possible estimates of the continuum energy levels, new simultaneous fits of the two-point 
functions in the Hg, Gig and G2g irreps were performed, in which the fitted energies for the matching ^ = | and 
J = I levels were forced to be equal: 



J^(Hg) 


77i(G2g) 
— ^1 ' 




77i(G2g) 


EiHg) 
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Continuum 



Splitting 



a ^ 0.11 Im 0.08 fm 



5 + 
2 




- si''"''' 


5.8(2.0) 


2.5(2.0) 


5 + 
2 




- 4^'' 


0.70(44) 


0.44(64) 


7 + 
2 






2.1(1.1) 


1.6(1.4) 


r+ 

2 




_^{Gl.,) 


1.49(78) 


0.38(79) 


7 + 

2 


p(G2g) 

-^2 


- 


0.59(45) 


1.24(72) 



TABLE V. Splitting of continuum energy levels with J > | into different irreducible representations of the double-cover 
octahedral group. All results in MeV. The data at a ~ 0.11 fm are from the ensemble with am„,d = 0.005, while the data at 
a ~ 0.08 fm are from the ensemble with amu,d ~ 0.004. 

This was implemented by augmenting the function of the simultaneous fit in the following way: 



E, 



(He, 



4 ~ ^1 



(G2g) 



E. 



^(G2.,) 



E, 



(He, 



/a^+ E 



-E. 



(G2g) 



2 



(33) 



where the energies E^i^^ are expressed in terms of the actual fit parameters as E^i^^ = bJ^-* 



<5(^^ + ... + ^ri- The 

width a in Eq. |33] was chosen about two orders of magnitude smaller than the typical statistical unc erta inty in the 
energies. By minimizing the augmented fit parameters are returned that satisfy the conditions (32 1 up to the 
input width a. These new fits still had x^/d.o.f w 1, because of the smallness of the energy splittings between the 
different irreps. Performing the simultaneous fit with the enforced relations (32 1 also stabilizes the extraction of the 



very close energy levels (such as i?g " and Ei^ ' ), and makes the spectral overlap factors more sharply peaked 



can be seen in Fig. [sj Note that in this work no further constraints beyond that of Eq. ( 33 1 were imposed on any of 
the fit parameters. 

These simultaneous fits, along with simple one-exponential fits in the and Gi„ irreps, yield 11 different bbb 
energy levels. Having performed the angular momentum identification, these levels can now be labeled by and a 
new subscript counting the states in each channel by increasing energy: 





E2{\^ 






E2{f 






E2ir 










EiiD, 






EiiD- 







(34) 

Because NRQCD is used in this work, the extracted energies do not include the rest masses of the three b quarks, 
i.e. they are all shifted by a common amount that is not known with sufRcient precision. Therefore, only energy 
differences are considered in the following. 

The remaining point to be discussed in this section is the choice of tmin, the starting time slice from which the fits 
are performed. This parameter has to be chosen large enough such that the contamination from higher-excited states, 
which decay exponentially with t, is negligible. However, tmin must not be made too large either, as the statistical 
uncertainties increase with t^in and the fits eventually become unstable. Figures |6] and [7] show the tmin-dependence 
of the set of ten independent energy splittings chosen here. For the matrix two-point functions in the Hg, Gig, and 
G2g irreps, the total number of time slices included in the fit, Nt — tmax/a " imin/a + 1, was held constant as tmin 
was varied, to keep the dimension of the data correlation matrix fixed at a manageable size (Nt = 5, 8, 8 for the Hg, 
Gig, G2g irreps, respectively). 

As can be seen in Figs. [g] and [jj for the energy splittings aEi (| ) — aEi ( § ^ ) , aE2 ( | ^ ) — o.Ei ( § ^ ) , and aEi ( 5 ^ ) — 

ai?i(|^), which are large energy differences between bbb states of rather different spatial structure, the plateaus set 
in later than for the other, smaller splittings, which mainly constitute the fine- and hyperfine structure. The values 
of ^min that were chosen to get the best estimates of the energy splittings are given in the captions of Figs. |6] and [7] 
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FIG. 6. Dependence of the results for the bbh energy splittings on the start time slice tmin that is used in the fit. The data 
shown here are for the ensembles with a ~ 0.11 fm, with the light quark masses of am^.d ~ 0.005, 0.01, 0.02, 0.03 (from left to 
right). The shaded bands indicate the best possible estimates of the energies, which are taken from tmin/a = 8 or tmin/a. = 7 
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all other energy splittings. 
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FIG. 7. Dependence of the results for the bbb energy splittings on the start time slice tmin that is used in the fit. The data 
shown here are for the ensembles with a « 0.08 fm, with the light quark masses of amu,d ~ 0.004, 0.006, 0.008 (from left to 
right). The shaded bands indicate the best possible estimates of the energies, which are taken from tmin/a = 12 for the three 
large energy splittings a_Ei(| ) — aEi{^'^), aiJ2(|^) — ai?i(|'*'), aiJi(|'^) — a_Ei(|^), and from tram /a = 6 for all other energy 
splittings. 
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V. FINAL RESULTS FOR THE bbb SPECTRUM 



In the previous section, ten bbb energy splittings were computed for each of the seven different ensembles of gauge 
fields. These results are given by the horizontal bands in Figs. |6] and [7] The values of the light sea-quark masses used 
in the generation of the gauge field ensembles correspond to pion masses that are larger than physical (see Table IV I . 



The final step of the analysis is to perform extrapolations of the bbb spectrum to the physical value of the pion mass. 
These extrapolations are done here using the same method that was used for the bottomonium spectrum in Ref. [37] . 
The light quarks influence the bbb spectrum only through their vacuum-polarization effects, and the dependence on 
mu,d is weak. Therefore, it is sufficient to perform the extrapolations linearly in rriu.d, and hence linearly in m^. 

The bbb energy splittings were first converted to MeV using the values of the lattice spacings as given in Table 
JV) Then, coupled fits to the data for the two different values of the gauge couphng, /3i = 2.25 and /32 — 2.13, were 
performed using 



(35) 



where E{in'^,/3) denotes a generic bbb energy splitting. The ensembles with /3 — f3i have a « 0.08 fm, while the 
ensembles with P = f32 have a ~ 0.11 fm. The free fit parameters in Eq. (35) are E{0,Pi), E{0,f32), and A. No 



continuum extrapolation is performed here, because lattice NRQCD is an effective field theory that requires a cut-off 
^ mil- The only assumption made here is that higher-order effects proportional to terms like a?rn^ are negligible, 
so that the same parameter A can be used for both values of (3. 

The fits to the data for the ten bbb energy splittings using Eq. (351 are visualized in Fig. [Sj Evaluating the fitted 
functions for m^r = 138 MeV leads to the results given in Table [Vll In addition to the ten independent energy splittings 
discussed so far, the Table also gives some further combinations for convenience, in particular the energy differences 
of all ten excited states to the ground state i?i(|^), and a result for the very small splitting £'4(|^) — i?2(|^) that, 
as a consequence of the strong correlations, has smaller absolute uncertainties than the other splittings involving the 
same levels. 

As can be seen in Fig. [8] and Table |VI[ the results for the bbb spectrum show only a weak dependence on the 
lattice spacing, which in most cases is not statistically significant. The results at a « 0.08 fm and — 138 MeV 
can be quoted as the predicted values for the continuum bbb spectrum, once the remaining systematic uncertainties 
have been estimated. These estimates can be made using information from Sec. |VI| about the dependence of the bbb 
energy splittings on the couplings Ci in the NRQCD action [see Eq. (24)]. The systematic uncertainty is computed 
individually for each energy splitting E, using the formula 



(syst) _ 



dE\^ 



da) 



(0.02^si)' + (0.07 {E-Esi) 



-I 1/2 



(36) 



which takes into account the varying contributions from spin-dependent and spin-independent NRQCD interactions. 

The first two terms in Eq. ( 36 ) correspond to the uncertainty in E that results from the uncertainty in the tuning 
of the NRQCD coefficients C3, and C4 [see Eq. 



27[)]. The derivatives with respect to C3 and C4 are approximated 
using discrete difference quotients formed from the results in the last three columns of Table |VII| To save computer 
time, the results in Table [Vll] were obtained at the coarser lattice spacing a ~ 0.11 fm. However, for the purpose of 
estimating a^^^^^\ it is sufficient to approximate the derivatives with respect to C3 and C4 at a 0.08 fm as being 
equal to those at a ■ 



0.11 fm, and then setting CTcs = 0.084 and (Tc4 = 0.053 according to Eq. (27 1 for a « 0.08 fm 



The third term in Eq. (36) describes the systematic uncertainty in the spin- independent contribution to the energy 
splitting. This contribution, Egi, is obtained by setting C3 = C4 = C7 = cs = cg = in the NRQCD action. Given the 



weak a-dependence of the spectrum, Egi can be taken from the second column of Table VII However, the estimate of 
a 2% systematic uncertainty is specific to a « 0.08 fm. It includes the radiative, discretization, and relativistic errors, 
and is based on the discussion of radial and orbital energy splittings for the same lattice spacing in bottomonium |37j . 
The estimates of uncertainties for bottomonium are also valid also for triply-bottom baryons, since the energy- and 



momentum scales involved are the same (indeed, the results of Sec. VI confirm that the ?; -expansion converges at a 



similar rate for the bbb system as for bottomonium). 

The last term in Eq. (36) describes the systematic uncertainty in the spin- dependent contribution to the energy 
splitting. This contribution can be isolated by computing the difference {E — Esi), where E is the result from the full 
NRQCD action. Because the leading spin-dependent couplings C3 and C4 have been tuned nonperturbatively (and their 
tuning uncertainty is already taken into account), and because the spin-dependent order-w^ terms have been included 
in the NRQCD action at tree-level, the dominant remaining sources of error for the spin splittings are discretization 
errors and the missing radiative corrections in the u^-terms. Following the discussion of the bottomonium fine- and 
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TABLE VI. Energy splittings in MeV between various bbb states, extrapolated to the physical pion mass. In the final results 
(last column), the central values and statistical uncertainties are taken from a ~ 0.08 fm, and estimates of the total systematic 
uncertainties computed using Eq. ( 36 1 are given. The ground-state mass is equal to -Bi(f + ) = 14371 + 4+ 11 MeV [25]. 



hyperfine splittings in Ref. [37] , a systematic uncertainty of 7% is assigned here to the spin-dependent contributions at 
a « 0.08 fm. Again, the values of {E — Esi) can be taken from Table VII (the differences of the results from columns 
six and two), because the spectrum has a weak a-dependence. 

The final results for the bbb spectrum, with systematic uncertainties computed using Eq. (36 1, are given in the last 
column of Table VI The energy differences of the ten excited states to the ground state ^tbb are plotted in Fig. |9| 
The results for the different energy levels are highly correlated, and the small splittings between nearby states can 
in fact be computed with much smaller absolute uncertainties. These smaller energy splittings are given in the lower 



part of Table VI and are plotted in Fig. 10 



It is interesting to compare the QCD results obtained here to the potential-model calculation of Ref. [8] (see Fig. 5 
therein). The number of states in the considered energy region are in agreement, and the energy differences to the 
ground state predicted by Ref. [8^ are found to be within 10% of the QCD results. However, the potentials used 
in Ref. [S] did not include any spin-orbit or tensor interactions, so that the results obtained there have the exact 
degeneracies £^2(5 ) 



£;3(| + ) = Ei{^^) = Eiil'^), ^4(1"^) = E2{^^), and Ei^ ) = Ei{^ ). As can be seen in 



Fig. 10 the QCD calculation performed here is so precise that the spin-dependent effects that lift these degeneracies 



VI 



are clearly resolved. These effects will be discussed further in Sec. 

Reference [5^ also calculated the higher-lying bbb spectrum, and these additional states were all found to be separated 
by energy gaps of order 300 MeV from the states considered here. Along with the plateaus observed in Figs. [6] and [Tj 
the large energy gaps found in Ref. [8 provides further confidence that the contamination from higher states in the 
fits of Sec. |IV| is negligible. 

Remarkably, the three energy splittings £^2(5^) — -^1(5 ) — and £'2(|^) — Ei{^^) that were 

computed in the early bag-model calculation of Ref. [3] also agree with the results obtained here to within 10%. On 
the other hand, the energy splittings calculated recently using a quark-model in Ref. fTB] (see Table 19 therein) are 
in dramatic disagreement with the QCD results obtained here: by about a factor of two for the larger splittings and 
by about a factor of 10 for the smaller splittings. 
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FIG. 8. Extrapolation of the bbb energy splittings to the physical pion mass. The fits are linear in ml, and were done 
simultaneously for the data at the two different lattice spacings. The data are plotted with closed symbols, and the extrapolated 
results at = 138 MeV are plotted with open symbols. The fitted functions and their 1-sigma statistical uncertainty are 
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FIG. 9. Final results for the bbb spectrum relative to ground state -Ei(| 

values). The superimposed shaded regions show the statistical/scale setting and the total uncertainties, respectively, 
results are highly correlated, and the uncertainties for energy differences between nearby states are in fact much smaller than 
suggested by this plot. See Fig. 10 for close-ups of the spectra near and -Bi(| ), where advantage of the correlations 



is taken by computing the energy differences relative these levels 
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in the vicinity of these levels. The superimposed shaded regions show the statistical/scale setting and the total uncertainties, 
respectively. See the last column of Table \VJ\ for the numerical values. 
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VI. DEPENDENCE OF THE SPECTRUM ON THE COEFFICIENTS IN THE NRQCD ACTION 



In Sec. |Vj the bbb spectrum was computed with coefficients Ci in the lattice NRQCD action tuned such that the 
effective field theory reproduces relativistic QCD. Table [Vl| and Figs. [9l[l0| give the best possible results obtained here 
for the bbb energy levels in the real world. However, with lattice NRQCD, one can perform simulations for arbitrary 
values of the coefficients The ability to selectively turn on and off the different terms in the NRQCD action and 
compute the effect on the bbb energy levels can be exploited to gain deeper insight into the interactions between three 
heavy quarks. 

The numerical results of this section are summarized in Table |VII[ Shown there are the values of the bbb energy 
splittings computed for eight different choices of the coefficients in the NRQCD action. The various terms in the 
NRQCD action were already discussed in Sec. Ill and their coefficients Ci were defined in Eq. (24). The calculations 
in this section were done for a single gauge field ensemble only (a w 0.11 fm, am^^d = 0.005), to save computer time. 
As shown in Sec. |Vj the dependence of the bbb spectrum on a and rriu^d is weak, and therefore a single ensemble is 
sufficient for the purpose of studying the Ci-dependence. In all cases, the b quark mass and the Symanzik-improvement 
coefficients in the NRQCD action remained unchanged {ami, — 2.487, C5 = cg = 1). The following discussion focuses 
on the energy regions near ) and ), as this is where all the spin-dependent level splittings are found. 



The energy splittings in the first column of Table VII were computed with the order-w^ NRQCD action, which 
contains only Hq = — ^^A^^-' (and the associated lattice discretization improvement terms with C5 and cq). Turning 

on also the spin- independent order-w'' terms, — cig^ (A^^^)^ and C2g^ (V ■ E — E ■ V), gives the results in the 

second column of Table |VII| These results are plotted in Fig. [TT] In both cases, the action does not depend on the 
heavy-quark spin, so that L and S become separately conserved quantum numbers, up to the small effects of rotational 
symmetry breaking introduced by the lattice. In the absence of rotational symmetry breaking, one would then have 

the exact level degeneracies £^2(5"^) = £^3(1^) = £i(|^) = £1(5"^), £^4(1^) = £^2(1"^), and £1(5") = £i(i"). The 

relations Ei{^ ) = £'i(| ) and £2(5^) — ^sll^) actually remain exact on the lattice, an observation that can be 



related to the trivial subduction of these two J values into lattice irreps (cf. Sec. II B). The degeneracies with J > 



are only approximate, but the splittings remain very small. Note that the energies quoted here for the higher- J levels 
were obtained by averaging over the different irreps into which a continuum level splits [see the discussion around 
Eq. ( 33 ) ; also see Table |V] for the size of the original splittings between the different irreps] . 

Next, Figure 12 shows the spectrum after additionally turning on the leading interaction with the chromomagnetic 
moment of the heavy quark: 



C4 CT • B. 

2mi, 



(37) 



subtr 



This interaction causes small positive splittings [-£^2(^ ) — )] 

2.23(66) MeV, [£i(i^) - £i(i^)],,btr. = 2-05(52) MeV, and [^4(1^) 
rotational-symmetry-breaking-induced splittings seen at C4 = (second column of Table VII I have been subtracted 



: 1.45(84) MeV, [£;3(- 

^2(|^)],ubt 



2 ' '^1(2 )] subtr. 

= 4.28(47) MeV, where the 



The operator (37) also introduces a very significant splitting of the two odd-parity levels considered here: i?i(| ) 



) = —12.97(34) MeV. For heavy quarkonium, the operator iSln is mainly associated with spin-spin and tensor 
interactions. However, simple potential models for baryons that include only spin- spin and tensor interactions predict 
Ei{^ ) — £i(| ) = [52H54] . Thus, one can conclude that the operator (37) also plays an important role in the 



generation of spin-orbit interactions. This can indeed be seen in the derivation of spin-dependent potentials using 
pNRQCD [55 . 

The other spin-dependent interaction of order is given by 



- C3 



8ml 



(^ 



X E - E 



(38) 



Setting C4 = again, and turning on the interaction ( |38[ ) instead, produces the results shown in Fig. |13| For the 
bbb levels consi dere c 

introduced by (37): 



bbb levels consi dere d here, the operator (38) results in spin splittings with the opposite sign compared to those 

-18.63(42) MeV, [^3(1^) 

r5 + ^ 



ibtr 



-15.58(36) MeV, 



[£1(1 - £;i(r)]subtr. = -8.89(42) MeV, Ei[if) - £2(1 )]subtr. = -8.74(24) MeV, and Ei{^ 



i-^-£4(i- 



5* = |, the effect of 



7.05(13) MeV. Notice in particular that for the bbb levels with approximate structure L = 2 
(38) is an order of magnitude larger than the effect of (37). Furthermore, the shifts introduced for these levels by the 
operator (38 ) are approximately proportional to 2 i • S' = J{J -I- 1) — L{L -|- 1) — S{S + 1). This is what is expected 
for a spin-orbit interaction in baryons levels with totally symmetric spatial wavefunctions |56j . 
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TABLE VII. Dependence of the bbb spectrum on the coefficients Ci in the NRQCD action [see Eq. (24l]. All results are given 
in MeV. The data are from the ensemble with a ~ 0.11 fm and arriu d = 0.005. 
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Next, Fig. Il4| shows the bbb spectrum with both (jSTl) and ([38|) turned on (fifth column of Table Ivill). For [£^2(5"^) - 



turning on ( 37 1 and ( |38[ ) 
the splitting [-E'i(f "^^ 



and Eij^ ) - -Eld )' new results are consistent with the sums of the results from separately 



but there is some evidence for nonlinear behavior in the other spin splittings. For example, 
2 ) — -Bi(^~'^)] is equal to —4.17(59) MeV now, while the sum of the splittings obtained from 
separately activating (37) and (38 1 is —6.85(67) MeV. Of course there is no reason to expect linearity here: the lattice 



calculation is fully nonperturbative. 



Having included both (37) and (38), the action is now complete through order v . As can be seen by comparing 



the results in the first and the fifth columns of Table VII the radial and orbital bbb energy splittings obtained with 
the order-u^ and ordei-v'^ NRQCD actions differ by ^ 10%, demonstrating the convergence of the NRQCD expansion 
0.1 as in bottomonium. Finally, turning on additiona lly t he spin-dependent order-« ^ te rms by setting 

6 



with V « 

cr = cs = Cg — 1 gives the results in the sixth column of Table |VII[ which are plotted in Fig. [15 



The order- 



terms affect some of the bbb spin splittings by as much as 30%, showing that including these terms is essential to 
obtain precise results. Most of bbb spin splittings considered here decrease in magnitude when the order-w^ terms are 
included in the NRQCD action, as is familiar from bottomonium [37]. However, one notable exception to this rule is 
found here: the order-w^ corrections increase the mag nitude of i;i(i )-£^i(| )• 
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FIG. 11. Dependence of the spectrum near E\{1^^) and E\{\ ) on the coefficients Ci in the NRQCD action (at a ~ 0.11 fm, 
o.mu,d ~ 0.005). Shown here is the case of the spin-independent order-u* NRQCD action, obtained by setting cs = C4 — cr = 
cs ~ cg — 0. In the absence of rotational symmetry breaking, this leads to the exact degeneracies -B2(|^) ~ -Ead^) = = 



+ E4{^~^) = ^2(1^), and ) = ). On the lattice, the relations Ei{^ ) = EiC^ ) and £2(|^) = EaCf 



still exact, but the degeneracies with J > f levels are only approximate. 
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FIG. 12. Dependence of the spectrum near _Ei(|^) and ) on the coefficients Ci in the NRQCD action (at a ~ 0.11 

fm, arriu.d ~ 0.005). Shown here is the case of the order-i;* NRQCD action, but with the coefficient of the operator 

(T ■ ( V X E — E X V ] set to zero, so that the only remaining spin-dependent interaction is — C4 cr ■ B. 
\ / 2mb 
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FIG. 13. Dependence of the spectrum near and ) on the coefficients Ci in the NRQCD action (at a ~ 0.11 fm, 

amu,d = 0.005). Shown here is the case of the order-w* NRQCD action, but with the coefficient of the operator er ■ B set to 
zero, so that the only remaining spin-dependent interaction is 
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FIG. 14. Dependence of the spectrum near ) and ) on the coefficients d in the NRQCD action (at a ~ 0.11 fm, 

a-rnu,d ~ 0.005). Shown here is the case of the complete order-u'' NRQCD action. 
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FIG. 15. Dependence of the spectrum near 51(1 + ) and E-id ) on the coefficients Ci in the NRQCD action (at a ~ 0.11 fm, 
o.mu,d = 0.005). Shown here is the case of the complete NRQCD action as used in the main calculations of this work, including 
all terms of order-w* as well as the spin-dependent order-u® terms. 



VII. CONCLUSIONS 



In this work, the first nonperturbative QCD calculation of the baryonic analogue of the bottomonium spectrum 
was performed. By combining improved lattice NRQCD [SS] with other powerful techniques that have been developed 
more recently, the energies of ten bbb excited states were computed with high precision. The calculations include 2+1 
dynamical flavors of light quarks, and the bbb spectrum was extrapolated to the physical pion mass. The main results 
are given in Table [Vl] and are plotted in Figs. [9| and [TOl 

The reliable identification of triply-bottom baryon states with angular momentum up to J = ^ was greatly simplified 
by using interpolating operators constructed with the subduction method of Ref . [34 . As already observed in Ref. [34] 
for light baryons, the cross-correlations between interpolating operators subduced from different values of J are small. 
In the present work, it was additionally shown that these overlaps decrease when the lattice spacing is reduced. 
Furthermore, it was possible to resolve the small energy splittings of continuum bbb levels with J > | into the 
different irreducible representations of the double-cover octahedral group. It was shown that these splittings also 
decrease when the lattice spacing is reduced (see Table |v| , providing another demonstration of rotational symmetry 
restoration. While the suppression of mixing between different J-values is a general consequence of the approximate 
rotational symmetry, additional suppressions were observed here for the triply-heavy baryon two-point functions 
between operators constructed using different values of L or S. This feature is likely to be a consequence of the large 
b quark mass, resulting in a suppression of the spin-orbit coupling and hence an approximate individual conservation 
of L and S (the total orbital angular momentum and total quark spin). 

To implement the b quarks on the lattice, an NRQCD action including the spin-dependent order-w^ terms was 
used here, and the coefficients of the spin-dependent order-w"* terms were tuned nonperturbatively. Together with 
the high statistics, this allowed the calculation of the bbb spin splittings with ~1 MeV total uncertainty. To learn 
more about the forces between three heavy quarks, additional simulations were performed on one ensemble for several 
"unphysical" choices of coefficients in the NRQCD action, thereby disentangling the contributions of different NRQCD 
operators to the bbb energy splittings. These additional simulations also clearly demonstrated the convergence of the 



velocity expansion for bbb baryons, and facilitated the estimates of the systematic uncertainties given in Table VI 



The lattice QCD results obtained here for the triply-bottom baryon spectrum provide a unique opportunity to test 
quark models for baryons in the regime were the description using potentials is expected to work best. Most of the past 
potential-model calculations of baryon excited states have focused on light baryons, for which some experimental data 
are available. However, quark-model descriptions are bound to remain poor approximations for these complicated 
systems. Now that precise lattice QCD results for the much cleaner bbb spectrum can serve as a substitute for 
experimental data, it is desirable to perform new continuum-based calculations for triply-heavy baryons, using for 
example the quark model of Ref. [57j , or the modern pNRQCD approach [211 [13] . 
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